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ABSTRACT 



This research was conducted to enhance understanding of the use of 
high explosives to simulate the effects of a nuclear underwater explosion. A 
review of the known characteristics of the nuclear, spherical conventional, 
and tapered conventional underwater pressure-time histories illustrates the 
selection of the tapered charge to simulate the underwater nuclear explosion. 
Three areas of study were then pursued. The first compared the structural 
response resulting from attack by conventional and nuclear type pressure 
profiles, verifying the need to match duration as weU as peak pressure when 
simulating the underwater nuclear explosion. The second employed finite 
element analysis to study the three dimensional shock generated by a 
tapered charge. Third, a computer program was written to couple an 
optimizer with an existing tapered charge pressure-profile generating code to 
improve the tapered charge design process. 
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I. INTRODUCTION 



This thesis addresses efforts to improve simulation, using conventional 
high explosives and scale models, of the underwater shock environment and 
structural effects resulting from the underwater nuclear explosion. 
Simulation using scale models and small conventional charges provides 
valuable information without the requirement for nuclear testing, with 
minimal environmental impact, and at low cost. Better understanding of the 
physics involved impacts ship and weapons designs. 

In the current atmosphere of reduced military spending in order to reap 
the benefits of the "Peace Dividend" resulting from the end of the cold war, 
the threat to U. S. Navy ships and submarines from underwater nuclear 
explosion would appear to be greatly reduced. However, two factors ensure 
the continued existence of the threat from underwater nuclear explosion: 

1. The presence of the former Soviet Union's vast arsenal of nuclear 
weapons combines with economic instability to increase the Likelihood of 
more nations gaining access to the material necessary to construct 
nuclear weapons. 

2. The ceaseless march of technology’ worldwide dictates future 
growth in the number of nations attaining the particular technology 
necessary to build and detonate nuclear weapons. 

The growing community of nuclear weapons capable nations may not posses 

the same restraint from the use of nuclear weapons displayed since 1945. 

Given the existence of a threat from underwater nuclear explosion, the 
significance of this threat can be determined only if the effects are well 
understood. This same understanding is essential to incorporation of shock 
hardening in ship and submarine designs to improve survivability. 
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With the overall objective of improving ship and submarine 
survivabihty through better understanding of the phenomena associated 
with underv\mter explosions, the Naval Postgraduate School conducts 
ongoing research into underwater explosions and effects. This thesis is the 
result of a part of that continuing research, the first at the school related 
specifically to the nuclear underwater explosion. 

The known characteristics of the nuclear underwater explosion, 
together with a discussion of modeling techniques is found in Chapter II. 

Chapter III presents a comparison of the structural response of a simple 
cylinder subjected to side-on attack by conventional spherical charge and 
nuclear type pressure profiles. The doubly asymptotic boundary assumption 
combined with an explicit finite element method was used to perform the 
analysis. Results verify the need to match peak pressure and duration of the 
nuclear pressure-time history when designing test charges to simulate the 
underwater nuclear explosion. 

The tapered charge, due to inherent long pressure duration, commonly 
generates the simulated nuclear pressure field used in model testing. 
Chapter IV explores the three dimensional aspects of the shock front 
generated by a conventional underwater tapered charge explosion. The 
results of analysis performed using an exphcit finite element method are 
presented with the intent of supplementing existing knowledge of the 
tapered charge pressure profile which to date consists mostly of on-axis data. 
The results show the evolving shape of the shock front at early times and 
provide information on the relationship between the peak pressures 
measured on and off the charge axis. 



2 



As outlined in Chapter V, a computer code was written to optimize the 
design process for a conventional tapered high explosive charge. Starting 
with a desired pressure-time history and initial estimates for charge 
geometry and standoff distance, this program utilizes public domain 
optimization software to return improved design values. Although tested 
with a subroutine based upon an existing routine for calculating the 
pressure-time history of a tapered charge, this program may be coupled with 
other existing or future routines for calculating the tapered charge 
pressure-time history. 

Conclusions and recommendations for further research in this area may 
be found in Chapter VI. 
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II. CHARACTERISTICS AND MODELING OF THE 
UNDERWATER NUCLEAR EXPLOSION 



Before attempting to simulate the underwater nuclear explosion using 
conventional high explosive charges for scale model testing, a degree of 
familiarity with the shock generated by nuclear and spherical high explosive 
charges is warranted. This chapter, therefore, outHnes the use of empirical 
relationships to determine the pressure-time histories generated by the two 
types of charges. The tremendous weight and standoff distance required for 
a conventional spherical charge to create a pressure profile similar to that of 
a nuclear charge points directly to the need for scaling. 

Following a brief description of the principles used to construct scale 
models for underwater shock testing, these principles are apphed to nuclear 
and conventional spherical charges. The resultant large size and standoff, 
even after scaling, of the conventional spherical charge, leads to the selection 
of the tapered charge, made of conventional high explosives, to simulate the 
underwater shock generated by the underwater nuclear explosion. 

A. THE UNDERWATER NUCLEAR SHOCK PROFILE 

The energy content, or "yield", of a nuclear explosion is commonly 
measured in tons of TNT equivalent, the amount of explosive energy 
contained in 2,000 pounds of the conventional high explosive TNT. A one 
kiloton (kT) device contains the energy equivalent of 1,000 tons of TNT, one 
megaton (MT) that of 1,000,000 tons. Although the majority of the energy 
released in an underw^ater nuclear explosion contributes to the generation of 
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the underwater shock wave, the extremely high temperatures (tens of 
millions of degrees) reached in a nuclear explosion contribute to a significant 
amount of energy release in the form of thermal radiation. Chemical 
explosions, by contrast, occur at much lower temperatures (thousands of 
degrees), resulting in a higher percentage of the total energy released as 
kinetic energy to generate the underwater shock. (Glasstone and Dolan, 
1977, pp. 1-3,6, 11) 

To date, the United States has conducted fi.ve announced underwater 
nuclear explosions, from 1946 to 1962 (Bolt, 1976, pp. 251-274). Glasstone 
and Dolan (1977, pp. 268-272) provide three empirical charts to calculate the 
pressure-time history of an underwater nuclear explosion given the yield of 
the device and the standoff distance R from the explosion. From the curves 
of the first two charts, the maximum pressure (psi) and the time 
constant 0 (ms) are determined. The time constant equals, as in exponential 
decay, the time between the arrival of the shock when P = P^^^ and the time 
at which P = P^ax^^ ~ ^-^^Pmax- pressure actuaUy decays at a somewhat 
slower than exponential rate after one time constant, necessitating the third 
chart which plots the non-dimensional values P(t)/P„^,. vs t/6 to provide an 
idealized pressure-time history for an underwater nuclear explosion with no 
bottom or surface reflection efi'ects. Based upon these curves. Figure 1 shows 
the pressure time history of a 40 kT nuclear explosion at a standofi* of 1,000 
yds. Figure 2 illustrates standoff. Two prominent features of the 
underwater nuclear explosion stand out in Figure 1; the high pressure at a 
■sig nifi cant distance fi:om the explosion, and the long decay time. P^^^ 
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increases with yield, decreases with standoff: the duration of the shock wave 
increases with yield and standoff (Glasstone and Dolan, 1977, p. 269). 




Figure 1. Pressure profile of a 40 kT nuclear charge, E = 1000 yds. 



^ R 

G 

Charge 



Figure 2. Measurement of standoff distance R, 




Target 



In water of sufficient depth where bottom reflections are negligible or 
occur at a much later time, the steady pressure decay ends abruptly by the 
phenomenon of surface cutoff (Shin and Geers, 1991, §3.3; or Glasstone and 
Dolan, 1977, pp. 244-246). Figure 3 shows two paths followed by a shock 
wave emanating from an underwater explosion. The direct, compression 
wave strikes the target first with a sudden rise in pressure followed by the 
steady pressure decay described above. The other path shows a compression 
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wave striking the air- water interface and being reflected as a tension, or 
rarefaction wave. 



Air 




Approximating shock speed by the speed of sound in water gives the 
cutoff time t as the difference in distance traveled bv the direct and 

CO 

rarefaction waves divided by C^; 




The arrival of the rarefaction wave at time t^.^ after the arrival of the 
direct shock wave causes a sudden pressure drop as the remaining pressure 
from the compression wave is essentially canceled. Hence, if the 40 kT 



nuclear charge discussed earlier is detonated at a 285 ft depth in deep water, 
a target at 1,000 yards and the same depth will experience t^^ = 10.7 ms using 
the expression on the previous page with = 5 ft/ms. The resultant 
pressure profile with pressure cutoff is shown in Figure 4. 




Figure 4. Pressure profile of a 40 kT nuclear charge with surface cutoff. 
1,000 yds, D=285 ft. 



For an actual underwater nuclear explosion, the smooth curve of Figure 
4 would be further modified by refraction of the shock wave due to salinity, 
currents, and temperature variation in the water media. Whereas Figure 4 
depicts an instantaneous pressure rise upon arrival of the compression wave 
at time 0, a finite rise time would be expected. Additionally, the abrupt 
surface cutoff shown assumes the rarefaction wave and the compression 
wave travel at the same speed of sound in water. In actuality, the 
rarefaction wave travelling in shocked media partially overtakes the 
compression wave rendering a less steep pressure drop at t^^. For explosions 
occurring in shallow water, bottom reflections and retransmissions further 
alter the shock profile. Additional shocks may occur at later times due to the 
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bubble pulse resulting from explosions at depths such that the gases of the 
explosion expand and collapse before venting at the surface. More than three 
bubble pulses are unlikely due to steam condensation. (Glasstone and Dolan, 
1977, pp. 56, 245, 246, 269) 

Having established the general characteristics of the pressure-time 
history of the underwater nuclear explosion, the next section uses empirical 
formulas to determine the feasibility of emulating this profile using 
conventional high explosives of spherical shape. 



B. SPHERICAL CONVENTIONAL CHARGE PROFILE 

As in the case of the nuclear charge, the conventional charge of 
spherical or near spherical shape gives rise to a sudden pressure increase 
followed by exponenticd pressure decay for one time constant 0 and somewhat 
slower decay thereafter. Using exponential pressure decay for an 
approximation, empirical studies provide the follo^^dng useful relationships 
for determining the pressure-time history developed by a conventional high 
explosive of spherical or near spherical shape when detonated underwater 
(Shin and Geers, 1991, §3.2); 



P(t) = Pm*.xe 



t 

0 



P 



max 



= K 



0 = KsW 




= maximum pressure (psi) 

A2 

= time (ms) for P to decay to 



P max 

e 



where Ki , K 2 , Ai , A 2 = empirical constants for a given explosive 

t = time (ms) 

P = pressure (psi) 

R = standoff (ft) 

W = charge weight (lb). 
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Using the empirical relationships above, a charge of 27.5 thousand tons 
of TNT would be required to emulate the pressure- time history generated by 
a 40 kT underwater nuclear explosion as illustrated in Figure 5. 




t (ms) 



Figure 5. Comparison of 40 kT nuclear (solid line) and 27,500 ton TNT 
(dashed line) explosions. R=3,000 ft, D=285 ft for both charges. 



The fact that 27.5 kT of TNT generates a similar pressure profile as a 
nuclear device of 40 kT yield is attributed to the thermal energ>’ released by 
the nuclear explosion as discussed in Section II.A. Obviously, using 27.5 
thousand tons of TNT to study the effects of a nuclear explosion is not 
practical, nor feasible. This amount of high explosive represents a volume 
greater than that of the Washington Monument (Lapp, 1980). Scaling laws, 
discussed in the next section, enable the use of scale models and charges of a 
more reasonable size to simulate the underwater nuclear explosion and its 
effects. 
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C. HOPKINSON SCALING 



Dimensional analysis yields scaling principles which enable the use of 
scale models to repHcate the behavior of fuU size objects, or protot.vpes. Of 
particular utiUty in the analysis of underwater explosions and effects is 
Hopkinson scaling. Through the use a scale factor X, the quantities length, 



mass, and time are scaled as follows; 
Model 


Prototype 


XL 


Length 


L 


Xi 


Time 


t 


X^M 

giving the invariant quantities: 


Mass 


M 


P 


Density 


P 


a, P 


Stress, Pressure 


o, P 


8 


Strain 


8 . 



Some quantities, chiefly the hydrostatic loading due to gravity, are not 
adaptable to scaling and require additional consideration in modeling. (Shin 
and Geers, 1991, §4.1) 

The benefits of scaling are many. By reducing standoff and charge size, 
tests can be conducted in small manmade ponds wdth little or no 
environmental impact. Geometrically similar models can be constructed of 
the same materials as the prototype to study structural response to 
underwater shock. Model stress and strain levels ^\dll match those of the 
prototype. Since the material required to build model charges and structures 
is equal to X^ that of the prototype, model testing delivers obvious cost 
benefits. 
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Using a scale factor of = 1/30 to simulate the 40 kT nuclear 
pressure-time history discussed previously would require a spherical TNT 
charge of size 

Wm = X'^Wp = ^ X 28.5 x 10^ x tons x2, 000^ = 2, 040 Ibm. 

V 30 V ton 



The standoff required is 



Rm = >.Rp = 000 ft = 100 ft. 

oU 



Figure 6 shows the 40 kT nuclear pressure-time history scaled using X = 1/30 
together with the pressure profile of 2,040 pounds of TNT. The only 
difference between this figure and Figure 5 is the time scale. 




t (ms) 



Figure 6. 2,040 Ibm TNT charge (dashed line), R=100 ft, to simulate a 
40 kT nuclear (solid line) explosion using a scale factor of X = 1/30. 



Although the spherical conventional charge size is now feasible, the 
weight and standoff required for the simulation are too great for small pond 
testing limited to nominal charge weights and standoffs in the tens of pounds 
and feet respectively. One might be tempted in scale model testing to use a 
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smaller charge and shorter standoff to match only the peak pressure of the 
nuclear profile when studvdng shock effects on a scale model. As shown in 
Figure 7, a 14 pound spherical TNT charge at a standoff of 19 ft gives the 
same maximum pressure as a 2040 pound charge at a standoff of 100 ft. The 
pressure of the lighter charge decays much more rapidly due to its smaller 
time constant. 




Fi^re 7. Pressure-time histories of a 2,040 Ibrn (solid line) TNT charge 
at R=100 ft with a 14 Ibm (dashed line) TNT charge at R=19 ft. 



It \\dU be shown in Chapter 111 that the structural response from two 
pressure profiles having the same peak pressure but different duration of 
high pressure give rise to dramatically different structural responses. 
Therefore, in order to simulate the slow pressure decay of the underwater 
nuclear explosion using conventional charges of modest size, the shape of the 
conventional charge must be modified. The resulting shape is that of the 
tapered charge discussed in the next section. 
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D. TAPED CHARGE PRESSURE PROFILE 



Designed to generate long duration shock waves to simulate nuclear 
underwater shock loading on scale models, tapered charges consist of a series 
of truncated cones on a common axis fitted with a detonator on the small or 
nose end as shown in Figure 8. Constructed in sizes ranging from a few 
ounces to over 15,000 pounds, the tapered charge generates a directional 
pressure field wdth maximum duration along the nose side on the charge 
axis. (Gordon and Davidson, 1983) 




For tapered and spherical charges of the same weight and at the same 
standoff. Figure 9 illustrates the trade-off between peak pressure and 
duration distinguishing the two designs. The tapered charge creates a peak 
pressure at a lower value than the spherical charge, followed by a more 
gradual pressure decline over a region ctdled the pressure plateau. At the 
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end of the pressure plateau, the pressure-time history of the tapered charge 
gives way to exponential decay then surface cutoff. 




Figure 9. Pressure profiles of tapered (solid line) and spherical (dashed 
line) charges with the same weight and standoff. 



The pressure plateau generated by the tapered charge gives an obvious 
advantage over the spherical charge in simulating the nuclear shock profile. 
The duration of the pressure plateau t^, can be determined to a first 
approximation by estimating the time difference between (1) travel from nose 
to tail along the charge axis a distance L at the detonation speed Cp then a 
distance L -t- R to the target at the speed of sound in water Cg , and (2) travel 
from the nose to the target a distance R at speed Cg . The resulting 
expression for pressure plateau duration is 



tp - 



L L + R 
Cd Cs 



_ ^ - T ^ J_ ^ J-^ 
Cs ^Cd Cs^ 
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The above formula predicts a plateau duration independent of standoff. As 
Figure 10 shows, however, the duration of the pressure plateau actually 
decreases somewhat with increasing standoff. 




Figure 10 also shows the approximate 1/R relationship between peak 
pressure and standoff. In addition to varying the standoff to match peak 
pressures and the overall length to match pressure plateau duration, the 
segment lengths and diameters of the tapered charge can be varied to taylor 
the shape of the tapered charge pressure-time history to model a particular 
nuclear pressure profile. Design of a tapered charge involves first examining 
existing data for a charge design which produces a pressure-time history 
most nearly matching that desired. The charge design is then adjusted using 
computer calculations tempered wdth the experience of the designer. Small 
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scale experiments may be used to verify the design before production of a 

I 

I large tapered charge. (Costanzo, 1991) 

[ Additional knowledge of the tapered explosion as well as any 

I 

streamlining of the design process should improve the design of tapered 
charges used to simulate the underwater nuclear explosion. An apphcation 
of the finite element method (FEM) to study the early time propagation of the 
shock generated by a tapered charge explosion is found in Chapter IV. 
Chapter V outlines the apphcation of computer design optimization 
techniques to enhance the tapered charge design process. 
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III. EFFECT OF PRESSURE DURATION ON 
STRUCTURAL RESPONSE 



As mentioned in Section II. C, if one could match only the peak pressure 
when conducting a simulation of the underwater nuclear attack, a spherical 
charge of modest size would suffice. As illustrated in this chapter, however, 
matching of peak pressures alone is not adequate. Comparison of the 
response of a simple cylindrical shell to side-on attack by a long duration, 
nuclear t>q)e, and a short duration, conventional type pressure profQe yielded 
substantially different results using numerical techniques. 

A. ATTACK CUR\T:S 

In order for a realistic comparison in the model testing environment, the 
pressure profiles used in this study were generated from two 56 Ibm HBX-1 
charges. The short duration pressure profile was derived from the empirical 
relationships of Section II. B for a spherical charge at a standoff of 20 ft. The 
long duration pressure profile used in this study to simulate the nuclear 
profile was derived from scaled tapered charge data. Figure 11 shows the 
two attack curves \Nith the same peak pressure of approximately 3,400 psi 
but very different high pressure duration times. Each pressure-time history 
mth corresponding standoff was entered into an existing computer code to 
provide underwater shock loading of a simple cylindrical aluminum shell 
model as described in the next section. 
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3500 




Figure 11. Two attack curves used for comparison of structural 
response. 



B. NUMERICAL MODEL 

The attack geometry for this study is shown in Figure 12. The figure 
illustrates the standoff for a spherical charge. Tapered charge standoff is, by 
convention, measured from the nose, Figure 8, to the target. The pressure 
profiles and standoffs from the previous section were used \\dth a side-on 
attack cylindrical shell model developed by Fox (1992. pp. 60,61). The 
analysis was conducted using the public domain finite element method 
(FEM) code VEC/DYNA3D (Hallquist and Stillman, June, 1990) coupled \^dth 
the boundary element method code USA (Deruntz, 1989). The USA 
(Underwater Shock Analyzer) code reduces the media surrounding the 
cylinder and the associated forces to discrete forces and masses to provide 
loading to the cylinder. The FEM code VEC/DYNA3D utilizes exphcit time 
integration to provide the response of the cylinder to the apphed shock 
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loading. The coupling of the two codes was initiated by the Naval 
Postgraduate School (Fox, Kwon, and Shin, 1992). 





Charge 



Fi^re 12. Attack geometry for cylindrical shell FEM experiment. 



The target consisted of a 1/4" thick 6061-T6 aluminum cylinder with 1" 
machined and welded endplates. Material properties for the cylinder were: 

Oy = 40,000 psi 

E = 10,800 kpsi 
v = 0.33 
p = 174 Ibm/ft^ 

Figure 13 depicts the 550 shell elements comprising the one-quarter 
symmetry discretization of the cylindrical shell model used in this study. A 
description of the theory behind the shell elements used can be found in the 
article by Belytschko, Lin, and Tsay (1984). Symmetric boundary conditions 
were applied on the yz and zx planes. The boundary element loading 
described previously pro\dded the boundarj" conditions for the outer surface 
of the wet elements of the aluminum cylinder. 
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Figure 13. Quarter symmetry FEM model of simple cylindrical shell. 



The aluminum was modeled as a kinematic/isotropic/elastic/plastic. In 
this idealization, elastic deformation occurs at stress levels lower than the 
yield stress mth plastic deformation occurring at the yield stress. This 
provides a reasonable approximation of the uniaxial stress-strain curve for 
6061-T6 aluminum based upon pull-test data as shown in Meyers and Murr 
(1981, p. 40), with the exception of no provision for failure once a maximum 
engineering strain of 7 to 9 per cent has been sustained. Taking 8 per cent 
plastic strain as a failure criterion and applying a factor of safety of 2, failure 
of the aluminum shell model was predicted for effective plastic strain in 
excess of 4 per cent Effective plastic strain is defined as (Ugural and 
Fenster, 1987): 



8p = 
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1/2 
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where e,, e^, and s., are the true plastic strain components. Since the 
maximum strain-rate calculated in this study was approximately 100 
in/in-sec, well below the 2000 in/in-sec required for appreciable strain-rate 
effects, (Meyers and Murr. 1981, p. 50), strain-rate hardening was not 
included in the model. 

C. RESULTS 

Figures 14 and 15 show the damage sustained by the cylindrical target 
as a result of attack by the nuclear and conventional pressure profiles with 
the same peak pressures. The quarter-symmetry model has been reflected 
about the yz and zx planes using the post-processor TAURUS (Hallquist, 
1990) to provide visualization of the full model in each case. Displacements 
for the conventional attack of Figure 14 were scaled by a factor of 10 to better 
display the damage pattern. No scaling was done in Figure 15 where the 
damage was much more extensive due to the nuclear type attack. 



element 214 
element 211 




EFFECTIVE PLASTIC 
STRAIN 



2.798E-03 
5.595E-03 
8.393E-03 
1. 120E-02 
1 . 400E-02 



Figure 14. Effective plastic strain from conventional attack (displace- 
ments scaled by a factor of 10). 
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element 214 
element 211 




EFFECTIVE PLASTIC 
STRAIN 

2. 000E-02 
3. 250E-02 
4. 500E-02 
5.750E-02 
7. 000E-02 



Figure 15. Effective plastic strain from nuclear attack. (No scaling of 
displacements) 



Figure 16 illustrates the time histor>' of eflfective plastic strain for the 
elements sustaining maximum damage from each attack. The maximum 
damage of 1.67 per cent effective plastic strain experienced by the 
conventional attack, although significant, did not exceed the failure criteria 
selected in the previous section. The effective plastic strain from the nuclear 
type attack, exceeding 4 per cent after less than 0.6 ms, indicates the 
prediction of catastrophic failure of the cylinder resulting from this attack 
based upon the numerical analysis. 

Not only is the magnitude of effective plastic strain much greater for the 
nuclear type attack, the mode of damage is quite different. In both cases, 
maximum damage occurred in elements located approximately 3.2 inches 
from the endplate. In the conventional case, the maximum damage occurred 
in element 211, on the yz plane. Maximum damage in the nuclear case, 
however, occurred at element 214, 30 degrees off the yz plane. 
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however, occurred at element 214, 30 degrees off the yz plane. The locations 
of the two critical elements are shown on Figures 13 through 15. 



C 

’(0 

u. 

C/) 




Figure 16 . Plot of effective plastic strain for the elements sustaining 
maximum damage. 



Based upon the analysis of this chapter, both the peak pressure and 
high pressure duration must be generated by the conventional high explosive 
used to simulate the underwater nuclear explosion. Due to the ability of the 
tapered charge to accomplish this task at a lower charge weight than the 
spherical charge, the studies of Chapters IV and V were conducted to 
enhance understanding and improve the tapered charge design process. 
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IV. THE TAPERED CHARGE SHOCK FRONT 



Due to its ability to develop a long duration pressure profile from a 
modest charge weight, the tapered charge commonly provides the shock 
loading to simulate the underwater nuclear attack as discussed previously. 
In order to better understand the shock developed by tapered charges, finite 
element analysis (FEA) was conducted to study the shock developed by a 
tapered charge detonated underwater. The FEA provided information on the 
directional nature of the pressure-time history generated by the tapered 
charge. This information, specific for the charge geometry and t>T)e of 
explosive, can be used to determine the accuracy required to position the 
charge and target used in model studies of the underwater nuclear explosion. 

A. FEA MODEL 

Figure 17 shows the 46" x 46" x 152" quarter-symmetry FEA model used 
for this study. The mesh consists of 120 HBX-1 charge elements and 40.692 
water elements. Appendix A provides a Listing of the input file for the 
INGRID (Stillman and HaUquist, 1991) mesh generating program used. The 
dimensions of the tapered charge, corresponding to those of Figure 8, were: 

Li =L2 = 0.333 ft 
U = 4.333 ft 
L = 5.000 ft 
di = 1.125 in 
d 2 = 2.625 in 
ds = 4.125 in 
d 4 = 5.375 in. 
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Using a cast specific gravity for HBX-1 of 1.71 (Dobratz, 1981, p. 19.53) gives 
the charge weight of 60.1 Ibm. Referring to Figure 17, symmetric boundary 
conditions were applied on the yz and zx planes. Non-reflective boundary 
conditions were applied to the remaining four planes. 




Fi^re 17. FEA model used to examine the three dimensional aspects of 
the tapered charge. 



The charge and water model were analyzed using FEA program 
VEC/DYNA3D (Hallquist and Stillman. 1990). This is the stand-alone 
version of the coupled program used in Chapter III. The explosive was 
modeled as a high explosive bum material with the Jones-Wilks-Lee (JWL) 
equation of state, the water as a null material using the Gruneisen equation 
of state. 
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For the charge, the high explosive burn material model used by 
VEC/DYNA3D requires entry of the detonation velocity D, the Chapman- 
Jouget pressure P^j, and the density p (Hallquist and Stillman, 1990, p. 40). 
The values D = 0.731 cm/ps, P^j = 0.2204 Mbar, and p = 1.712 gm/cm^ for 
HBX-1 were taken from Dobratz (1981, p. 19.53). The JWL equation of state 
was used to describe the pressure-volume-energy behavior of the detonation 
products (Dobratz, 1981, p. 8.21); 

Where A, B, and C = linear coefficients in Mbar 
R] , R 2 , and co = nonlinear coefficients 
_ V _ volume of detonation products 

volume of undetonated high explosive 
P = pressure in Mbar 

E = detonation energy' per unit volume in • 

cm'^ 

The parameters E, A, B, Rj, R^, and co listed above, empiricaUy derived from 
cyhnder-test data, are required entries for use of the JWL equation of state 
with VEC/DYNA3D (HaRquist and Stillman, 1990, p. 89). For this study 
these values were taken from cylinder-test data for H-6, an explosive of 
similar composition to HBX-1, due to non-availabihty of HBX-1 data. The 
values used were (Dobratz, 1981, p. 8.22): 

E = 0.103^^^212^ 

cm" 

A = 7.5807 Mbar 
B = 0.08513 Mbar 
R] =4.9 
R 2 = 1.1 

CO = 0.20. 
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For the water, the VEC/DYNA3D null material model (Hallquist and 
Stillman, 1990, p. 41) requires entry of the density and an optional pressure 
cutoff. The value p = 1.000 gm/cm^ was used for density, and a pressure 
cutoff value of 6.89 x 10'® Mbar (0.1 psi) were used since water is unable to 
sustain tension. The Gruneisen equation of state wdth cubic shock velocity 
-particle velocity (u^-Up) defines pressure p in Khar for compressed materials 
as (Hallquist and Stillman, 1990, p. 91); 



P 



PoC^p 



1 + 



l-flM- 



Im- 





where 



p = density in 
Po 



gni 
cm3 

standard density in 



gni 

cm3 



and the required parameters are 

C = the intercept of the Us-u,, curve in 

Si , S2, S.3= coefficients of the slope of the u^-Up curve 
Yo - the Gruneisen gamma 
a = first order volume correction to yo . 

E = internal energy per unit volume in 



The parameters C, Sp S„, and S3, for the u^-Up curve were obtained from 
Steinberg (1987); 

U,=C + SiUp+S2^Up+S3(^)'up 

= 0. 148 + 2.56up - 1.986^Up + 0.2268[^J “up 
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Appendix B outlines the determination of = 0.4934 and a = 1.3937, based 
upon the seventh order polynomial approximation for the Gruneisen gamma 
as a function of specific volume developed by Gurtman, Kirsch, and Hastings 
(1971). 

B. RESULTS 

The large number of elements led to a computationally intensive finite 
element analysis of the tapered charge underwater explosion. Running on a 
UNIX engineering workstation required approximately four days to perform 
the computations through 2 ps after detonation. Figure 18 shows 
representative time histories of two elements located at a standoff of 35" from 
the charge nose. One element is on the charge axis off the nose, the other 90 
degrees off this axis. Due to the large variation in magnitude of pressures at 
the two locations, as will be shown in Section IV.B.2. the pressures have been 
normalized to better compare the general shapes of the curves. 

The element located on the charge axis maintained a greater portion of 
its maximum peak pressure over a longer period of time than did the element 
located 90 degrees off the axis. This observation further supports the use as 
well as the orientation of the tapered charge for simulation of the underwater 
nuclear explosion. 

Both of the curves of Figure 18 rose less rapidly than expected and 
displayed oscillatory behavior. Gordon and Davidson (1983) experienced 
similar results when analyzing a tapered pentohte charge using 
finite-difference techniques in two dimensions. They attributed the long rise 
time to three possible causes: the artificial viscosity coefficient built into 
their model, the mesh size, and the fact that their model was approximated 
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by a pointed nose. Although the model used in this study did not use a 
pointed nose approximation, the other two possible causes apply to this 
study. Fox (1992. pp. 9.10) cited mesh reflection effects as a contributing 
factor to instabilities when using the FEA code of this study. 




Figure 18. Normalized Pressure-time History of two elements located at 
R = 45 inches: on the charge axis and 90 degrees off the charge axis. 



Although desirable, a finer mesh was beyond the capabilities of the 
machine used in this analysis. For this study, the default artificial viscosity 
coefficient and time integration steps (Hallquist and Stillman, 1990, p. 25) 
were used. Future sensitivity studies may identify values for these 
parameters. The mesh reflection effect appears to be a strong factor in the 
observed pressure oscillations. The key contributor to this problem, uneven 
size of adjacent elements, was aggravated in this study by the necessity to 
conform the shape of eight-node elements to the curvatures of the tapered 
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charge. The end result was i&ne charge elements near the nose next to 
relatively coarse water elements and a reversal of this effect at the tail of the 
charge. 

Oscillations in pressure and the limited time of the pressure histories 
generated precluded determination of the pressure plateau duration. The 
remainder of this section concentrates on the shock wave as it emanates from 
the charge at early times and the directional nature of the peak pressures in 
the media surrounding the tapered charge. 

1. Early Time Shape of the Shock Front 

Figures 19 through 25 show a view perpendicular to the yz-plane 
of the water elements. The model has been reflected about the zx-plane, with 
charge elements removed. Figures 19 applies to time prior to detonation. 
Shading indicates the elements used for the plots of Section IV.B.2. 




Figure 19. Water elements prior to detonation. 
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Figures 20 through 25 comprise a set of pressure contour plots at 
early times through 0.58 ms. There are five contour lines in each plot. The 
contours range from 1,000 psi for the sparsest dotted line to 5,000 psi for the 
solid fine. This relatively low pressure remge at such close proximity to the 
charge was selected to clearly define the location of the evolving shock. 




Figure 20. Pressure contour plot of water elements at t — 0.08 ms after 
detonation. 



The effect of mesh reflection, discussed in the previous section, can 
be seen in Figure 20 along the charge axis to the front of the charge nose as 
the shock wave transmits through the media at different speeds in this 
region of very poor match of element sizes. The pressure contours bend 
inward toward the charge nose in the small element region. By 0.18 ms, 
Figure 21, the mesh reflection effect in front of the charge is negligible. 
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detonation. 
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detonation. 




Figure 24. Pressure contour plot of water elements at t — 0.48 ms after 
detonation. 



34 






Fig:ure 25. Pressure contour plot of water elements at t -0.58 ms after 
detonation. 



The shock front begins as a tear shape at early times emanating 
from the charge bum region and becoming nearly spherical on the charge 
axis off the nose. The above pressure contours show the nose portion of the 
evolving shock front to be nearly spherical in shape out to approximately 40 
degrees off the charge axis with a radius centered at the location of the 
undetonated charge nose. Once the bum region reaches the tail of the 
charge, the aft portion of the shock front expands to nearly spherical with a 
smaller radius of curvature than at the forward end. The above figures 
clearly illustrate the shock rapidly travelling outward ahead of the 
expanding detonation products. 
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2. Directionality of Peak Pressure 

This subsection examines variations in the peak pressure 
developed in regions surrounding the tapered charge based upon the FEA 
conducted. Figure 19 marks the locations of the elements whose 
pressure-time histories were gathered and processed to form the plots. 

For the plots of this section, relative peak pressure is defined 
as: 

n _ f* max 
x^rel - p 
^ o 

Where Pmax = maximum pressure computed for the element 

Po = maximum pressure observed for the element at the same 
standoff and nearest the charge nose axis. 

Figures accompanying the plots serve to further clarify the determination of 
P^p] as well as explain the geometry corresponding to each plot. 

Figure 26 corresponds to Figure 27, a plot of relative peak pressure 
as a function of 0 degrees off the charge nose axis from 0 = 0 to 90 degrees, at 
constant standofi', for three different standoffs. 




Figure 26. Constant standoff, 6 degrees off charge axis. 
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Figure 27. Relative peak pressure as a function of the angle off the 
charge axis at constant standoff. 



For the three standoff distances of Figure 27, there is a decrease in 
pressure from the on-axis pressure as the off-axis angle increases from 0 to 
20 degrees. This decrease, ranging from 3 per cent for R = 43 inches to 5 per 
cent for R = 15 inches, may be attributed to mesh geometry and mesh 
reflections. The effect of mesh geometry was caused by the centroid of each 
element not being at the exact standoff. This effect was minimized by 
applying a first order correction to each measured pressure: 

Ph = Ph4 

Where Pr = pressure used for plot at the nominal standoff R 
= max pressure for the given element at standoff R 

R = +y^+z'^ 

X, y, z = coordinates of the element centroid. 
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After approximately 20 degrees, the curves of Figure 27 rise 
steadily with increasing angle off the charge axis. This was to be expected 
due to the tradeoff between peak pressure duration and magnitude of the 
tapered charge. As 0 increases, the target is placed nearer to the bulk of the 
mass of the charge. In a test situation gages must often be located off the 
charge axis to record the free-field pressure experienced by the target 
without being overly influenced by proximity to the target. The test designer 
often attempts to locate the gage as far off axis as possible while minimizing 
the deviation of the gage measurement from the on-axis values. Figure 27 
indicate that, for the charge and standoffs of this study, gages could be 
located up to 35 degrees off the charge axis for an error in maximum peak 
pressure measurement of less than 5 per cent. To keep the error in 
measurement of plateau duration to an acceptable level, the allowable 
off-axis angle may lie well below the level based upon peak pressure alone. 

Past 90 degrees off the front charge axis, constant standoff was 
replaced by constant distance y off the charge axis to examine the relative 
pressures experienced by elements to the side of the charge. Figure 28 
illustrates the geometry involved. for this case is taken at a standoff R = y 
on the axis off the charge nose as shown in the figure. Figure 29 shows 
relative peak pressure as a function of ^ / L off the charge axis for ^ / L = 0 to 
1 from nose to tail. As can be seen from the figure, relative peak pressure 
shows marked departure for different values of y. While the highest relative 
peak pressure occurred at ^ / L = 0.8, the magnitude fell from 8.4 at y = 19 
inches to 6.2 at y = 43 inches. This was likely a close-in phenomenon. It is 
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Figure 28. Constant distance off the charge axis, to the side of the 
charge. 




Figure 29. Relative peak pressure as a function of the fraction of the 
length aft of the nose at constant distance off the charge axis to the side. 



expected that, at distances further removed from the charge axis than those 
of this study, the maximum relative peak pressure to the side of the charge 
will continue to fall then reach a steady value significantly lower than 
calculated here. Of interest is the fact that the maximum relative peak 
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pressure observed for elements located to the side of the charge was found at 
^ / L = 0.8, significantly aft of the center of gravity location at ^ / L = 0.6. 
This may be due to a build-up in pressure away from the center of gravity 
location in the direction of the tapered charge high explosive burn. 

Completing the examination of pressures encircling the charge 
from fore to aft, Figure 30 shows the geometry corresponding to Figure 31 
plotting the variation in relative peak pressure at a constant standoff R 
measured from the center of the tail. The angle 0 in this case is measured 
from 90 degrees off the rear charge axis toward the rear charge axis in order 
to continue proceeding in a counterclockwise direction. As in the two other 
cases, is measured off the charge nose, in this case at a standoff equal to 
the constant R. As in the constant R case off the charge nose, a correction 
was applied to the calculated values to account for variation in standoff of the 
element centroids. 




Figure 30. Constant R off the rear charge axis. 

Figure 31 shows a continuation of the decrease in peak pressure 
for the standoffs studied until 0 = 60 degrees, or 30 degrees off the charge tail 
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Fi^re 31. Relative peak pressure as a function of degrees aft of the 
charge. 



axis followed by a sHght increase to the charge axis. Here again, the 
unevenness of the curves near the axis may be attributed in part to mesh 
reflections. As in the side case, the relative pressure curves, though 
exhibiting similar shape, varried significantly for different standoffs, with 
shorter standoff corresponding to higher relative peak pressure. Though of 
minor interest in nuclear simmulations, the side and rear relative pressure 
plots are of more interest to weapons and industrial high explosive designers. 
The higher pressures at the tail end of the charge agree with the use of this 
configuration for demolition and other applications. 

Of more interest in simulations is the effect of distance off axis at a 
constant standoff from the charge nose as illustrated in Figure 32. This 
configuration, somewhat easier to set up for experimental verification, 
yielded the plot of Figure 33 which provides information of similar utility as 
the constant R plot of Figure 27. 
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Figure 32. Constant distance from nose along charge axis. 




Figure 33. Relative peak pressure at constant distance along charge 
axis as a function of degrees off the charge axis. 



Figure 33 indicates that, for the studied charge design and stand 
off distances, the off axis distance may be varied up to 10 degrees before 
exceeding a 5 per cent error in peak pressure readings. Since the allowable 
angle to stay \vdthin a given error tolerence decreased with range, a smaller 
angle would be allowed at the longer ranges used in simulations of the 
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underwater nuclear explosion. Having concluded the FEA of the underwater 
tapered charge explosion, the next chapter outlines computer optimization of 
a simple method to determine the pressure-time history of a tapered charge. 
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V. TAPERED CHARGE DESIGN OPTIMIZATION 



This chapter outlines the coupling of a public domain computer 
optimization package ADS (Vanderplaats, 1984) and a simple tapered charge 
pressure-time generating program TAPER (Costanzo, 1991) to optimize the 
tapered charge design process. The coupling program, listed in its entirety in 
Appendix C, can be used with existing and future tapered charge 
pressure-time history generating codes to enhance the design of the tapered 
charges used to simulate underwater nuclear explosions. 

In addition to the program in Appendix C, ADS (Automated Design 
Synthesis) and a pressure-time generating subroutine are required to 
perform tapered charge design optimization. 

A. SIMPLE PRESSURE-TIME HISTORY ALGORITHM 

Simpler, less computationally intensive, computer codes than that used 
for the finite element analysis of the previous chapter exist to predict the 
pressure profile of an underwater tapered charge explosion. These simpler 
codes are usually based upon an empirically based superposition principle. 

Because the empirically based exponential approximation of Chapter II 
works well for spherical / near spherical charges, one method of deriving the 
pressure-time history of an underwater tapered charge explosion is to 
partition the tapered charge into subsegments as shown in Figure 34. Each 
subsegment is considered to be a separate charge generating its own 
pressure-time history. (Costanzo, 1991) 
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Figure 34. Tapered charge discretized, into subsegments for calculation of 
pressure-time history using summation method. 



Shocks generated by subsegments detonated after the nose subsegment 
travel in shocked media at faster speeds than preceding shocks, overtaking 
them. There is thus a summation or stacking effect of individual "wavelets" 
to determine the overall pressure-time history of the tapered charge 
underwater explosion (Costanzo. 1991): 

A -Sl 

1=1 

Where P = pressure as a function of time 
N = number of waves stacked 
A, = empirical amplification factor, a function of 0i 
0i = slope angle of subsegment 
ri = R (^i = subsegment standoff 
R = standoff from charge nose to target 
^ = distance from nose to subsegment midpoint 
B = empirical decay constant 
t = time 

Di = charge diameter at subsegment midpoint 
5^ = subsegment length . 

Figure 35 illustrates the superposition scheme for a simple tapered charge 
discretized into large segments. 
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Figure 35. Siinimation of wavelets generated by charge subsegments to 
obtain resultant pressure-time history for a tapered charge. 



One algorithm, using a superposition scheme and empirical data to 
determine the pressure-time history of a tapered charge detonated 
underwater is the PASCAL program TAPER written by Fred Costanzo 
(1991). This program was converted to a FORTRAN program then to the 
FORTRAN subroutine used for tapered charge design optimization. Besides 
a basic overhaul of the input / output, a simple provision for surface cutout 
was added to the original algorithm. 

B. OPTIMIZATION PROBLEM 

The objective of the tapered charge design optimization was to 
determine charge geometry, as depicted in Figure 8, and standoff required to 
develop a pressure-time history most nearly matching a desired pressure- 
time history. The average square root of the sum of the squares of the 
differences between the computed and desired pressures at each time 
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increment was selected to quantify the difference between the optimal design 
and desired pressure profiles. Constraints restricted the charge to increasing 
diameters and placed upper and lower limits on charge weight, charge 
dimensions, and standoff. Only designs within these limits were allowed. 
For a tapered charge of three segments, the optimization problem became: 



k(p,-P,) 

minimize: J S 77 , i = 1, N 

\ i=l 

subject to: 1 in < L < 20 ft, i = 1 to 3 

0.75 in < dj < 10 in, i = 1 to 4 

5 ft < R < 30 ft, 

25 Ibm < W < 125 Ibm. 
dj < d, < d3 < d, 

where N = number of computed pressure-time history points 

Pj = desired pressure at a particular time 
Pj = computed pressure at a particular time 
= charge segment length 
d^ = charge joint diameter 
R = standoff 
W = charge weight. 

The main FORTRAN program DTAPOPT was written to accomplish the 
tapered charge design optimization using, as described previously, the ADS 
optimization package and the subroutine based upon the TAPER program. 

C. OPTIMIZATION PROGRAM 

In addition to DTAPOPT, three subroutines for use in conjunction with 
the main program were also written. The first, DCOMPAR, computes the 
square root of the average sum of the squares of the pressure differences. 
The second, DPTGEN, interpolates to find the desired pressure 
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corresponding to the design pressure at a given time. This subroutine 
enables the comparison of the desired and design pressures at the same 
times. The third subroutine, DTAPWT, computes the weight of the tapered 
charge for evaluation of the weight constraint. 

1. Required Program Input 

DTAPOPT requires input from two files and the keyboard. Figure 
36 shows a sample of a tapin.dat file used to input initial charge geometry, 
standofi", nominal length of pressure-time history to be computed, and surface 
cutoff time. Decimal alignment and horizontal placement on the lines is 
optional with one number per line only. Number translations have been 
written into the figure. 



3 


tmmber of charge segments 


0.5 


length of first segment (ft) 


1.0 


length of second segment (ft) 


5.0 


length of third segment (ft) 


1.0 


first joint, or nose, diameter (in) 


2.0 


second joint diameter (ifi) 


3.0 


third joint diameter (in) 


5.0 


fourth joint, tail, diameter (in) 


15.1 


standoff (ft) 


1.12 


nominal lerjgth of computed time history^ (ms) 


1.095 


surface cutoff time (ms) 



Figure 36. Sample tapin.dat file. Italics, not part of the file, indicate 
the meaning of each nu mber. 

Figure 37 shows a sample of the second input file, profin.dat, used 
to input the desired pressure profile for program DTAPOPT. Up to 1001 data 
points may be entered without program modification. The first time entered 
must be at time zero, and the last time must be greater than the nominal 
length of computed time history entered in tapin.dat. Decimal alignment 
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and horizontal placement on the hnes is optional. Only one number can be 
used on the first fine, two on the rest. Number translations added to the 
figure are shown in itahcs. 



9 iuIps, InUiJ iwml>pr of dodrod press} iro-timr data points niitws 1 



0.0 


0.0 


tdes(0)m ms-must be zero. 


pdes(O) hi psi 


0.04 


1950.0 


tdcs( 1) 


pdes(I) 


0.119 


1755.0 


tdes(2) 


pdes(2) 


0.278 


1560.0 


tdes(3) 


pdes(3) 


0.417 


1365.0 


ides(J) 


pdes(i) 


0.635 


1170.0 


ide.s(5) 


pde.sO)) 


0.834 


975.0 


ides(6) 


pdes(6) 


1.072 


780.0 


tdes(7) 


pdes(7) 


1.120 


0.0 


tdes(3) 


pdes(S) 


3.0 


0.0 


tdes(9) 


pdes(9) 



Ides(ndes) must Wf<reulrr than nominal length ofp-t history entered in lapin.dat 

figure 37. Portion of a sample profin.dat file. Italics, not part of the 
file, have been added to indicate what each number represents. 



Three integers comprise required keyboard input. These numbers 
control the optimization method used by ADS. Vanderplaats (1985) provides 
more detailed instruction on method selection. Vanderplaats (1984) provides 
the theory behind the methods. The three numbers consist of any one 
number from each of three groups. The first number may be any of the 
foUot\ing to determining the optimization strategy’'; 



First 

Number Optimization Strategy 



0 Go directly to the optimizer 

6 Sequential Linear Programming 

7 Method of Centers (Design must be feasible) 

8 Sequential Quadratic Programming 

9 Sequential Convex Programming. 
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The second required input number may be either of the following to 
determine the optimizer: 



Second 

Number 


Optimizer 


4 

5 


Method of Feasible Directions 
Modified Method of Feasible Directions. 



The last input number, one of the following, selects the one-dimensional 
search option to be used: 



Third 

Number 


One-Dimensional Search 


5 

6 

7 

8 


Golden Section Method 
Golden Section Method Plus 

Polynomial Interpolation 
Bounded Polynomial Interpolation 
Unbounded Polynomial Interpolation. 



2. Program Output 

Output from DTAPOPT includes one screen summarizing the 
optimization and an output file containing the design and interpolated 
desired pressure-time histories. 

Figure 38 shows a sample screen from a run of DTAPOPT on a 
personal computer. The figure shows the execution commands followed by 
prompts for the three inputs, in this case the user selected the combination 
8-5-7 for the strateg>% optimizer, and search. The output summary then lists 
the initial and final charge designs as well as the average square root of the 
pressure differences for each case and the number of calls to the pressure- 
time generating routine. 
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C : \P6cC\F0R\0PT : dtapopt 
ENTER ISTRAT.IOPT.IONED: 857 



INITIAL DESIGN 








LENGTHS = 


1.000 


1.000 


5.000 


DIAMETERS - 


1.00 


2.00 


4.00 


CHARGE WEIGHT = 


66.0 LB 






RANGE = 


20.0 FT 






S(?IT OF AVG OF (Pdesign-Pdesired)"2 


= 424.3 PSI 


FINAL DESIGN 








LENGTHS = 


0.083 


0.083 


4.825 


DIAMETERS = 


2.52 


3.75 


4.28 


CHARGE WEIGHT = 


55.7 LB 






RANGE = 


17.4 FT 






SQRT OF AVG OF (Pdesign-Pdesired)"2 


= 240.3 PSI 



FT 

5.00 IN 



FT 

4.52 IN 



CALLS TO P-T HISTORY GENERATOR = 158 



C:\P6cC\FOR\OPT: 

V 

Fi^re 38. Sample screen output from DTAPOPT. 




The output screen of Figure 38 represents an early step in the 
design process. The next would be to use this final design to input a more 
refined initial design, then run the optimization program again. For test 
optimizations, the square root of the average pressure difference squared was 
well below 100 for the "best" final design. 

Figure 39 shows a portion of a sample DTAPOPT output file. The 
first two columns contain the calculated time and pressure, the third column 
contains the interpolated values of the desired pressure profile input. The 
data are in free format to retain maximum precision, sacrificing the 
readability of formatting. This output file may be readily used to create plots 
comparing the initial and final designs as was done in the next section of this 
chapter. 
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0 . OOOOOOOOOOOOOOOOOE-01 


0 . OOOOOOOOOOOOOOOOOE-01 


0 . OOOOOOOOOOOOOOOOOE-01 


2 . 81022269002981981E-02 


1639 . 32749805998742000 


1405 . 11134501490983000 


5 . 620A4538005963962E-02 


1675 . 89437624031484000 


1987 . 35262142392480000 


0.47773785730506935 


1628.10858635031309000 


1658 . 35094063994575000 


0.50584008420536752 


1614 . 83610147863396000 


1636.41749525434739000 


0.53394231110566581 


1602 . 96393658242482000 


1614 . 48404986874857000 


1.06788462221133162 


900 . 78633^5776994000 


899 . 17511224684437800 


1.09598684911162980 0 


. OOOOOOOOOOOOOOOOOE-01 


235.12876911529235700 



Fi^re 39. Sample excerpted from an output data file generated by 



DTAPOPT. 



D. OPTIMIZATION RESULTS 

The design space using the subroutine adapted from program TAPER 
proved to be fraught with local minima attributable discontinuities resulting 
from integer changes in the number of waves stacked. Several runs of 
DTAPOPT using different optimization methods for a particular initial 
design resulted in an improved design unless the initial design was optimal. 
The improved design was then used as the initial design for further 
optimization. This process was repeated from four to seven times until 
further optimizations failed to improve the design an appreciable amount. 
Each run of DTAPOPT took approximately two minutes on a personal 
computer of modest, 386SX, capacity. Total time to perform a tapered charge 
design optimization was from one to two hours. 

Of the desired pressure profiles and initial designs tested, input 
parameter combinations 0-5-7, 8-5-7, and 9-5-7, usually produced the most 
improved design with the fewest calls to the pressure-time history generating 
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subroutine. These input numbers translate, as shown in Section V.C. 1 to 
direct optimization, sequential quadratic pro^amming, or sequential convex 
programming strategy" combined \vdth the modified method of feasible 
directions optimizer and a bounded polynomial interpolation 
one-dimensional search. For the best designs found, as mentioned in Section 
V.C. 2, the square root of the average difference between desired and design 
pressures was well below 100 psi. Figure 40 shows a comparison of the 
pressure profile generated by an optimized design from program DTAPOPT 
and the corresponding desired pressure profile. 




Figure 40. Pressure profile resulting from optimization compared with 
the corresponding desired pressure profile. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



The boundary element / finite element analysis of the side on attack of a 
simple cylinder by nuclear and conventional type pressure profiles resulted 
in substantially different structural responses. It is therefore necessary to 
model the duration of the pressure plateau to adequately simulate structural 
response to an underwater nuclear explosion. 

The finite element analysis of an underwater tapered charge explosion 
provided insight to the early time propagation of the shock wave and the 
directional variations in peak pressures developed in the surrounding media. 
There is need for future research, including sensitivity studies of the 
artificial viscosity coefficient and the time integration step, to obtain less 
oscillatory results to provide tapered charge pressure plateau duration 
information. Additionally, the computationally intensive method used lends 
itself more readily to supercomputer use where the mesh size may be 
extended far enough away from the charge form comparison with test data. 

Coupling of an optimization routine with a tapered charge pressure 
profile generating routine provides a tool which can be used to more 
efficiently design tapered charges used to simulate underwater nuclear 
explosions. The program written for this study may be used with existing 
and future pressure-time history codes. 
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APPENDIX A; TAPERED CHARGE FEM INPUT FILE 



Follo^^ang is a complete listing of the input file used to perform FEM 
analysis on the tapered charge and water model of Chapter IV: 



G0G02 

c Set up for DYNA3D, integrate to time 2000, TAURUS data dump interval 

10 . 

c high speed printer dump interval 29999. 

dn3d vec term 2000 plti 10.0 prti 29999.0 

c Scale X, y. z axes from inches to cm. 

xsca 2.54 ysca 2.54 zsca 2.54 

c Define material 1, type 8- -High explosive burn. HBX-1: 
c detonation vel 0.731 cm/us, Chapman- Jouget pres 0.2204 Mbar . 
c density 1.712 gm/cc . 

mat 1 type 8 d .731 pci .2204 ro 1.712 

c Equation of state for charge: 2--JWL: 

c a 7.5807. b .08513. rl 4.90. r2 1.10, omega .20, eO .103 Mbar-cc/cc. 

eos 2 a 7.5807 b .08513 rl 4.90 r2 1.10 omega 0.20 eO 0.103 endmat 
c Define material 2, type 9- -Null Material. Water: 

c cutout pressure .1 psi (6.89e-9 Mbar), density @ 4 deg C 1.000 gm/cc. 

mat 2 type 9 pc 6.89e-9 ro 1.000 
c Equation of state for water: 4- -Gruneisen : 

c sound speed .142 cm/us @ 4 deg C. si 2.56. s2 -1.986. s3 .2268, 
c Gruneisen gamma .4934, first order vol cor'n to gamma 1.3937. 

eos 4 sp .142 si 2.56 s2 -1.986 s3 0.2268 gamma .4934 sa 1.3937 endmat 

c Define two symmetry planes by defining a point and a normal vector 
c for each plane. Any point within .001 cm of a symmetry plane will be 

c included in that plane's definition. 

plane 2 

000 100 .001 symm 

000 010 .001 symm 

c Define detonation points in HBX-1. 

detp 1 point 000; 
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O 



c Define rear and front water cylinders, cone and plane surfaces for 
c charge nose (part 1) and transition water (part 2). 



sd 1 cn2p 00 0 001 0.5625 0.0 1.3125 -4.0 

sd 3 plan 00 0 001 

sd 4 plan 00-4 001 

c Part 1: small charge cone. 



start 
13 5; 

13 5; 

1 3; 

-10 1 
-10 1 
0.0 -4.0 

di 1 0 3; 103; ; 



sf i 


-1 


-3: 


; -1 


-3; 


: sd 1 






sf i 


-1 


-3: 


: -1 


-3; 


-1 -1; 


sd 


3 


sf i 


-1 


-3: 


; -1 


-3: 


-2 -2; 


sd 


4 


d 1 


1 1 


2 


3 2 










d 1 


1 1 


3 


2 2 











mate 1 
end 

c Part 2: water transition from small charge cone to square grid. 



start 



13579: 








13579; 








1 3; 








-4-4044 








-4-4044 








0.0 -4.0 








di 1 2 0 4 5 


: 1 2 0 4 5 ; ; 






d 2 2 0 4 4 


0 






sfi -2 -4 ; 


-2 -4 ; ; sd 1 






sfi -1 -5 ; 


-1 -5 : -1 -1 ; 


sd 


3 


sfi -1 -5 : 


-1 -5 : -2 -2 ; 


sd 


4 


dill 35 


2 






dill 53 


2 







mate 2 
end 

Define cone and plane surfs for med charge and transition (pts 3 

). 

sd 8 plan 0 0 -8 001 

sd 9 cn2p 00 0 001 1.3125 -4.0 2.0625 -8.0 

c Part 3: medium charge cone. 

start 
1 3 5; 

1 3 5; 

1 3; 

-10 1 
-10 1 
-4.0 -8.0 



and 
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di 1 0 3; 1 0 3: ; 
sfi -l-3:-l-3;:sd9 
sfi -1 -3: -1 -3; -1 -1; sd 4 
sfi -1 -3; -1 -3; -2 -2; sd 8 
dill 232 
dill 322 
mate 1 
end 

c Part 4: raediuni water transition. 

start 13579; 

13579; 

1 3; 

-4-4044 

-4-4044 



-4.0 -8. 


0 










di 1 2 0 


4 


5 


: 1 2 0 4 5 : ; 






d 2 2 0 


4 


4 


0 






sfi -2 


-4 




-2 -4 ; : sd 9 






sfi -1 


-5 




-1 -5 : -1 -1 : 


sd 


4 


sfi -1 


-5 




-1 -5 ; -2 -2 ; 


sd 


8 


dill 


3 


5 


2 






dill 


5 


3 


2 







mate 2 
end 

c Define cones and plane for large charge and transition (pts 5 and 6). 

sd 11 cn2p 00 0 001 2.0625 -8.0 2.6875 -60.0 

sd 13 plan 0 0 -60 001 

c Part 5: large charge cone. 

start 
13 5: 

13 5; 

1 27; 

-10 1 
-10 1 
-8.0 -60.0 



di 1 
sfi 


0 3:1 
-1 -3; 


0 3: : 
-1 -3; 


; sd 11 




sfi 


-1 -3; 


-1 -3; 


-1 -1; sd 


8 


sfi 


-1 -3; 


-1 -3; 


-2 -2; sd 


13 



dill 232 
dill 322 
mate 1 
end 

c Part 6: large water transition. 

start 

13579: 

13579 ; 

1 27; 

-4-4044 
-4-4044 
-8.0 -60.0 

di 12045:12045; ; 
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d 2 2 0 4 4 0 

sfi -2 -4 : -2 -4 : ; sd 11 

sfi -1 -5 ; -1 -5 ; -1 -1 ; sd 8 

sfi -1 -5 ; -1 -5 : -2 -2 ; sd 13 

dill 352 
dill 532 
mate 2 
end 

c Define cylinders for front and rear water cylinders and transitions. 

sd cl cyli 000 001 2.6875 

c Part 7: rear water cylinder. 

start 
1 3 5; 

13 5: 

1 24: 

-10 1 
-10 1 

-60.0 -106.0 
di 1 0 3; 1 0 3; ; 
sfi -1 -3: -1 -3: ; sdcl 
sfi -1 -3: -1 -3: -1 -1: sd 13 
dill 232 
dill 322 
nrb 002 002 

mate 2 
end 

c Part 8: rear water transition. 



start 

13579 ; 

13579 ; 

1 24: 

-4-4044 
-4-4044 
-60.0 -106.0 

di 12045: 12045; ; 
d 2 2 0 4 4 0 
sf 2 2 1 442 sd cl 

sfi -1 -5; -1 -5; -1 -1; sd 13 
d 111 352 

d 111 532 

nrb 0 0 2 002 

mate 2 
end 

c Define cones and plane for water cone and transition (parts 9 and 10) . 

sd 5 cn2p 000 001 0.5625 0.0 1.3125 4.0 

sd 7 plan 0 0 4 001 

c Part 9: front water cone. 



start 
13 5; 
13 5; 



58 



I 



I 



13; 

-10 1 
; -10 1 

\ 0.04.0 

) di 1 0 3; 1 0 3 ; ; 

' sfi -1 -3; -1 -3; ; sd 5 

I sfi -1 -3; -1 -3; -1 -1; sd 3 

1 sfi -1 -3; -1 -3; -2 -2; sd 7 

I dill 232 

j dill 322 

I mate 2 

end 



c Part 10: front water cone transition. 



start 
13 5 
13 5 
1 3; 

-4 -4 
-4 -4 
0.0 4 
di 1 2 0 4 
d 2 2 0 4 

sfi -2 -4 



9; 

9; 

4 

4 



4 

4 

5 

4 



sf i 
sfi 
d 1 
d 1 
mate 
end 



-1 -5 
-1 -5 
1 1 
1 1 
2 



; 1 2 
0 

-2 -4: 
-1 -5; 
-1 -5: 
5 2 
3 2 



0 4 5 



: sd 5 
-1 - 1 ; 
-2 - 2 : 



sd 3 
sd 7 



c Define cylinder for front water, parts 11 and 12. 
sd c2 cyli 000 001 1.3125 

c Part 11; front water cylinder. 



start 
13 5; 

13 5; 

1 22 ; 

-10 1 
-10 1 
4.0 46.0 

di 1 0 3; 1 0 3; ; 
sfi -1 -3; -1 -3; ; sd c2 
sfi -1 -3; -1 -3; -1 -1; sd 7 
dill 232 
dill 322 
nrb 002 002 

mate 2 
end 



c Part 12: front water cylinder transition. 



start 
13579 ; 
13579; 
1 22 ; 
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-4-4044 
-4-4044 

4.0 46.0 

dil2045: 12045 
<1 2 2 0 4 4 0 
sf 2 2 1 442 sd c2 

sfi -1 -5; -1 -5; -1 -1 
d 111 352 

d 111 532 

nrb 002 002 

mate 2 
end 

c Part 13: top water. 

start 
1 24; 

1 22 ; 

1 24 54 77: 

0.0 46.0 

4.0 46.0 

-106.0 -60.0 0.0 46.0 

nrb 2 0 0 200 

nrb 0 2 0 020 

nrb 0 0 1 001 

nrb 004 004 

mate 2 
end 

c Part 14: side water. 



start 
1 22 ; 

1 3; 

1 24 54 77: 

4.0 46.0 
0.0 4.0 

-106.0 -60.0 0.0 46.0 

nrb 2 0 0 200 

nrb 0 0 1 001 

nrb 004 004 

mate 2 
end 



end 



sd 7 
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APPENDIX B: GRUNEISEN GAMMA APPROXIMATION 



The seventh order approximation for the Gruneisen gamma as a 
function of the specific volume v developed by Gurtman, Kirsch. and 
Hastings (1971) is: 



y(v) = ao + aiv + a2V“ + • • -+a7v" 

where ao = 2, 366.6324 

ai =-22,669.420 
a2 = 91,259.368 
a.s= -200. 175.85 
ai =258,585.11 
a5 = -196. 872.84 
ae = 81,850.023 
ay =-14,342.530 

Y = Gruneisen gamma, dimensionless 

V = specific volume in ^^7. 



By definition. 



giving 






1 ^ 

po V 



V = 



Vo 

P+1 



P+1 



for p - 



g 

cm3 



X 

Vo* 



Substituting l/(p+l) for v into the Gurtman equation for p = 0 to 0.8 in 
.001 increments, then using a least squares linear fit with the same intercept 
resulted in the following hnear equation for y(p): 



y(p) = yo + ay = 0.4934 + 1.3937p 
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Fi^re 41 is a plot of y(li) from the Gurtman equation and the linear 



approximation. 



Grunciscn Gamma 




Fi^re 41. Comparison of y(p) from Gurtman, Kirsch, and Hastings 
equation with the linearized equation for y(p). 
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APPENDIX C: FORTRAN PROGRAM 



This appendix contains a complete listing of the FORTRAN main 
program DTAPOPT and its supporting subroutines DCOMPAR, DPTGEN, 
and DTAPWT, written to optimize tapered charge design. A separate 
subroutine to generate pressure-time histories given tapered charge 
geometry and standoff distance is required, as well as the optimization 
package ADS. 

PROGRAMMING NOTE; The program DTAPOPT and associated 
subroutines were written in double precision FORTRAN. The public domain 
version of ADS is written in single precision. ADS was converted to double 
precision for use with DTAPOPT. DTAPOPT and its three subroutines can 
be simply converted to single precision by removing the IMPLICIT NONE 
and DOUBLE PRECISION statements from the codes. 



I. PROGRAM DTAPOPT 



C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 



DTAPOPT 

This FORTRAN program combined with subroutines DCOMPAR, 
DPTGEN, and DTAPWT , is designed to optimize tapered charge design 
using the public domain optimization package ADS (converted to 
DOUBLE PRECISION by the author) 

--When coupled with a separate subroutine (not included) 
to calculate the pressure -time history of a three -segment tapered 
charge . 

Up to ten charge segments may be used with appropriate 
modifications to the input files. 



PRECISION: 

INPUT FILES: 

INTERRACTIVE INPUT: 
OUTPUT FILE: 



DOUBLE 

TAPIN.DAT Initial Design. 

PROFIN.DAT Desired P-T History. 
Optimization options. 

Named in TAPIN.DAT, plot data for computed 



time vs computed pressure and interpolated pressure. 

SCREEN OUTPUT: Starting and Optimized Designs. 

REQUIRED SUBROUTINES: DTAPWT Computes charge weight. 

DPTGEN Linearly interpolates to give desired 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



AUTHOR : 

MOST RECENT UPDATE: 
REFERENCES : 



pressure at each time output by P-T history generator. 

DCOMPAR Computes the square root of the 
average of the square of the difference between the computed and 
interpolated desired pressure at each time. 

ADS Package of subroutines which perform 

the optimization. 

P-T HISTORY GENERATOR. 

William Earl Miller II 
6/16/92 

(1) Vanderplaats , G. N., ADS - A FORTRAN 
Program for Automated Design Synthesis, Version 1.10. program 
instructions, Naval Postgraduate School, Montery, California, 

May. 1985. 

(2) Vanderplaats, G. N., Numerical 

Optimization Techniques for Engineering Design: With Applica- 

tions, McGraw-Hill Publishing Company, 1984. 

TEST P-T GENERATOR: The subroutine SDTAPER. was converted to 

FORTRAN by the author, based upon the PASCAL program TAPER version 
7/27/89 by F. A. Costanzo. 

DATA INPUT FILES- -Two required: 

TAPIN.DAT contains 12 lines, one value per line which are read 
into the program. The integer and decimal data need only be in 
a format suitable for list-directed input assignment to INTEGER 
and DOUBLE PRECISION data types respectively. The character data 
must be in proper form for a DOS file name. 



c 


LINE 




READ TO 


DATA 


c 


NUMBER 


DESCRIPTION VARIABLE 


TYPE 


c 


1 


number of charge segments 


NTAP 


integer 


c 


2 


length of charge segment (ft) 


TAP(l) 


decimal 


c 


3 


II 


TAP(2) 


decimal 


c 


4 


It 


TAP(3) 


decimal 


c 


5 


diameter of charge segment (in) 


DIAM(l) 


decimal 


c 


6 


If 


DIAM(2) 


decimal 


c 


7 


ft 


DIAM(3) 


decimal 


c 


8 


tt 


DIAM(4) 


decimal 


c 


9 


standoff (ft) 


RANGE 


decimal 


c 

c 


10 


nominal length of time history (ms^ 
used by P-T generating subroutine 


1 TLEN 


decimal 


c 

c 


11 


pressure cutout time (ms) 
used by P-T generating subroutine 


TCO 


decimal 


c 

c 


12 


name of output data file to be 
created 


NAME I L 


character 



c PROFIN.DAT contains NDES+2 lines, one integer value (NDES) on the 

c first line, two decimal values on each of the remaining lines, 

c Times must be in ascending order starting with 0.0 and extending 

c to a time greater than TLEN entered in TAPIN.DAT above. 



c 


LINE 




READ TO 


DATA 


c 


NUMBER 


DESCRIPTION 


VARIABLE 


TYPE 


c 

c 


1 


total number of desired pressure- NDES 

time history data pairs minus one 


integer 


c 

c 


2 first time first pressure 

(must be zero) 


TDES(O) PDES(O) 


decimal 


c 

c 


3 


second time second pressure 


TDES(l) PDES(l) 


decimal 


c 

c 

c 


NDES+2 


NDES time NDES pressure 


TDES (NDES ) PDES (NDES) 


decimal 



INTERRACTIVE KEYBOARD INPUT- -three integers required 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



These integers control the optimization method used by ADS. 

See Ref 1 for more complete instructions. Ref 2 for theory. 
Combinations 057, 857, and 957 are recommended. Others may 
work well in a given design space. Using the initial design, three 
or four runs with different combinations should result in a most 
improved design which can then be used as the initial design for 
further optimization. The three integers are: 



ENTRY 

ORDER 

first 



second 



third 



READ TO 
VARIABLE 
ISTRAT 



DATA 

TYPE 

integer 



DESCRIPTION 

optimization strategy used by ADS 
0 Go directly to the optimizer 

6 Sequential Linear Programming 

7 Method of Centers (Design must be feasible) 

8 Sequential Quadratic Programming 

9 Sequential Convex Programming 

optimizer to be used by ADS lOPT integer 

4 Method of Feasible Directions 

5 Modified Method of Feasible Directions 

one dimensional search options lONED integer 

5 Golden Section Method 

6 Golden Section Method Plus Polynomial Interpolation 

7 Bounded Polynomial Interpolation 

8 Unbounded Polynomial Interpolation 



OUTPUT FILE: Named in TAPIN.DAT. consists of NTIME+1 rows of 
three columns. Format is list-directed from DP variables. 

COLUMN WRITTEN FROM 

NUMBER DESCRIPTION VARIABLE 

one times output from P-T GENERATOR TIME(I), I=0.NTIME 

two P's from P-T GEN'R at each time PRESS(I). I=0,NTIME 

three int . P's for each time from DPTGEN PCOMPAR(I), I=0,NTIME 



SCREEN OUTPUT: Outputs initial and final design values plut the 

number of calls to the P-T History Generator. 

INITIAL AND FINAL VALUES 

DESCRIPTION VARIABLE 



lengths of tapered charge segments (ft) 

diameters of tapered charge (in) 

charge weight (lb) 

standoff (ft) 

sqrt of average sq diff of Pdesign-Pdesired 
OPTIMIZATION EFFICIENCY 



TAP(I) .I=1.NTAP 
DIAM(I) ,I=0.NTAP 
WEIGHT 
RANGE 
OBJ 



number of calls to P-T generator 



NCALLS 



VARIABLES 



NAME 


DATA 

TYPE 


COMMON 


ASSIGNED 

BY 


USED 

BY 


DESCRIPTION , 


ADS 

arg 

? 


A(I.J), 
1=1, NRA; , 


DP 

,J=1,NC0LA 


ADS 


ADS only 


Constr grads 


’ Y 


DF(21) 


DP 


- 


ADS 


ADS only 


Obj grads 


Y 


DIAM(I), DP 

I=1,NTAP+1 


OPTDPA 


TAPIN.DAT, MAIN 


MAIN. 

PTGEN'R 


Charge diaras 


N 


G(I), 


DP 


- 


MAIN 


ADS 


Constraints 


Y 


I=1,NC0N 


I 


I 


- 


MAIN 


MAIN 


Local indexing 


N 


IC(I), 


I 


- 


ADS 


ADS only 


Gradient ID 


Y 


I=1.NC0N 


IDG(I), 


I 


- 


MAIN 


ADS 


Constraint ID 


Y 


I=1,NC0N 
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c 


IGRAD 


I 


- 


MAIN 


ADS 


Grad Calc Ctrl 


Y 


c 


INFO 


I 


- 


MAIN. ADS 


MAIN, ADS 


Prog Flow Ctrl 


Y 


c 


lONED 


I 


- 


user 


ADS 


One-D Search 


Y 


c 


lOPT 


I 


- 


MAIN 


ADS 


Optimizer 


Y 


c 


IPRINT 


I 


- 


MAIN 


ADS 


ADS Print Opts 


Y 


c 


I STRAT 


I 


- 


user 


ADS 


Strategy 


Y 


c 

c 


IWK(I), 

I=1,NRIWK 


I 


“ 


ADS .MAIN 


ADS 


Work Array 


Y 


c 


NAMFIL 


CH 


- 


TAPIN.DAT 


MAIN 


Output File 
PTGEN'R calls 


N 


c 


NCALLS 


I 


- 


MAIN 


MAIN 


N 


c 


NGOLA 


I 


- 


MAIN 


ADS 


No . A columns 


Y 


c 


NCON 


I 


- 


MAIN 


ADS ] 


No. constraints 


Y 


c 


NDES 


I 


FORDPTGEN 


PROFIN.DAT 


DPTGEN Des^d plot sts-1 


N 


c 


NDV 


I 


- 


MAIN 


ADS ] 


No . Design vars 


Y 


c 


NGT 


I 


- 


ADS 


ADS 


Gradient Ctrl 


Y 


c 


NRA 


I 


- 


MAIN 


ADS 


A rows 


Y 


c 


NRIWK 


I 


- 


MAIN 


ADS 


IWK Dimension 


Y 


c 


NRWK 


I 


- 


MAIN 


ADS 


WK Dimension 


Y 


c 


NTAP 


I 


OPTI 


TAPIN.DAT 


PTGEN^R,MAIN 


No. Chg Segs 


N 


c 


NTIME 


I 


OPTI 


PTGEN'R DPTGEN . DCOMPAR Pit sts-1 


N 


c 


OBJ 


DP 


- 


MAIN 


MAIN ADS 


Obj Func Val 


Y 


c 

c 


PCOMPAR(I) , 
1=0, NDES 


DP 


FORDPTGEN 


DPTGEN 


DCOMPAR 


Interp'd Press 


N 


c 

c 


PDES(I) 
1=0, NDES 


DP 


FORDPTGEN 


PROFIN.DAT 


DPTGEN 


Desired Press 


N 


c 

c 


PRESS(I) 
1=0. NTIME 


DP 


OPTDPA 


PTGEN^R 


PCOMPAR,MAIN 


Calc Press 


N 


c 


RANGE 


DP 


OPTDP 


TAP IN. DAT. MAIN 


PTGEN'R 


Standoff 


N 


c 


SUMSO 


DP 


FORDCOMPAR 


DCOMPAR 


MAIN 


Press var'n 


N 


c 

c 


TAP(I) 
1=1, NTAP 


DP 


OPTDPA 


TAP IN. DAT, MAIN 


PTGEN'R, MAIN Chg Seg Lth 


N 


c 


TCO 


DP 


OPTDP 


TAPIN.DAT 


PTGEN'R 


Surf CO Time 


N 


c 

c 


TDES(I) . 
1=0. NDES 


DP 


FORDPTGEN 


PROFIN.DAT 


DPTGEN 


Des'd PT Time 


N 


c 

c 


TIME(I) 
1=0, NTIME 


DP 


OPTDPA 


PTGEN ^ R DPTGEN , DCOMPAR 


,MAIN Calc time 


N 


c 


TLEN 


DP 


OPTDP 


TAPIN.DAT 


PTGEN'R 


Lngth PT hist 


N 


c 

c 


VLB(I), 
1=1, NDV 


DP 


- 


MAIN 


ADS 


L lim on DV 


Y 


c 

c 


VUB(I), 
1=1, NDV 


DP 


- 


MAIN 


ADS 


U lim on DV 


Y 


c 


WEIGHT 


DP 


PASS 


DTAPWT 


MAIN 


Charge Weight 


N 


c 

c 


WK(I) . 
1=1, NRWK 


DP 


- 


ADS, MAIN 


ADS 


Work Array 


Y 


c 

c 

c 

c 


X(I), 

1=1, NDV 

SUBROUTINES 


DP 




MAIN, ADS 


ADS 


Des Vas 


Y 



c ADS ( INFO , I STRAT , lOPT , lONED , IPRINT , IGRAD , NDV , NCON , X , VLB , VUB , OBJ , 
c G,IDG,NGT, IC,DF,A,NRA,NCOLA,WK,NRWK, IWK,NRIWK) 

c See Ref 1 for more complete instructions, Ref 2 for theory, 
c Simply, ADS inputs design variables, constraints , and the objective 
c function for an initial design, then modifies that design, 
c requesting the corresponding objective and constraint values from 
c the calling program, 
c ARGUMENT VARIABLES 

c A(NRA.NGOLA) DP Array of constraint grads. ADS use only here, 

c DF(NDV+1) DP Array of objective gradients. ADS use only here, 

c G(NCON) DP Array of constraints for current design in X 

c G(l&2): 25#<WEIGHT<125#; G(3-5): DIAM(1)<DIAM(2)<DIAM(3)<DIAM(4) 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 



IC(NGT) 

IDG(NCON) 

I GRAD 

INFO 

lONED 

lOPT 

IPRINT 

ISTRAT 

IWK(NRIWK) 

NCOLA 

NCON 

NDV 

NGT 

NRA 

NRIWK 

NRWK 

OBJ 

VLB(NDV+1) 

VUB(NDV+1) 

WK 

X(NDV+1) 

X(I)=TAP(I) 

VLB(I)=1" 

VUB(I)=10' 

DCOMPAR 



I Array of constraint ID's. NA this program 
I Array ID'ing type of constraints: 0 for nonlin ineq 
I =0 for ADS calculate gradients using FD . 

I ADS flow control parameter, 

I =5 , 6 , 7 , 8 : See input section. 

I =4,5: See input section. 

I =0000 FOR NO ADS PRINTOUT 
I =0,6, 7, 8, 9: See input section 

I Stores ADS I vars . Some modifiable. 

I Dimensioned columns of A, min NDV+1 . 

I Number of constraints in G. 

I Number of design variables in X. 

I Returned to by ADS for gradients. 0 this program. 

I Dimensioned A rows: at least NDV+1. 

I Est : 200+NDV+NCON+N+MAX(N, 2^NDV) .N=MAX( NDV, NCOLA) 

I Es t : 500+10*(NDV+NCON)+NCOLA^^(NCOLA+3 )+N^(N/2 )+l 
DP Objective function value, SUMSQ provided by DCOMPAR 
DP Array of design variable lower bounds, indices as X 
DP Array of design variable upper bounds, indices as X 
DP Array for ADS double precision variables. 

DP Array of design vars. Assigned by input, then ADS. 
1=1,3; X(4+I)=DIAM(I) , 1=0,3 ; X(8)=RANGE 

VLB(4+I)=0.75'' VLB(8)=5' 

VUB(4+I)=10" VUB(8)=30' 



Communicates via COMMONS: OPTI , OPTDPA. FORDPTGEN , FORDCOMPAR 
Input variables: NTIME; PCOMPAR(I), PRESS(I), I=0.NTIME 

Output variable: SUMSQ 

DPTGEN 



Communicates via COMMONS: OPTI, OPTDPA, FORDPTGEN 
Input variables: NTIME; TIME(I), 1=0. NTIME 

NDES; TDES(I), PDES(I), I=0.NDES 
Output variables: PCOMPAR(I). 1=1. NTIME 
DTAPWT 

Communicates via COMMONS: OPTI, OPTDPA. PASS 

Input variables: NTAP ; DIAM( I ) , 1=0 , NTAP , TAP( I ) , 1=1 , NTAP 

Output variables: WIEGHT 

PRESSURE-TIME HISTORY GENERATOR 
Communicates via COMMONS: OPTI, OPTDP , OPTDPA 

Input variables: NTAP; RANGE; TLEN; TCO; TAP(I). 1=1. NTAP; DIAM(I). 

1=0 , NTAP . 

Output variables: NTIME; TIME(I), PRESS(I), 1=0, NTIME 



PROGRAM DTAPOPT 
IMPLICIT NONE 

SPECIFICATIONS FOR SUB VARIABLES 

COMMON/OPTI /NTAP , NTIME 
INTEGER NTAP, NTIME 
COMMON/OPTDP/ RANGE , TLEN , TCO 
DOUBLE PRECISION RANGE . TLEN , TCO 

COMMON/OPTDPA/ TAP, DIAM, TIME, PRESS 

DOUBLE PRECISION TAP(IO) , DIAM(0 : 10) ,TIME(0 : 1000) , PRESS(0 : 1000) 
COMMON/PASS/ WE I GHT 

DOUBLE PRECISION WEIGHT 

COMMON/FORDPTGEN/NDES , TDES , PDES , PCOMPAR 

INTEGER NDES 

DOUBLE PRECISION TDES(0:1000) , PDES (0 : 1000) , PCOMPAR( 0 : 1000) 

COMMON/FORDCOMPAR/SUMSQ 

DOUBLE PRECISION SUMSQ 
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on non non n on 



c 



c 



c 



c 



c 



SPECIFICATIONS FOR ADS VARIABLES 



up to 20 X's. 100 G's. 30 const grads 
INTEGER IWK(2000) ,IDG(100) . IC(30) 

DOUBLE PRECISION X(21 ) , VLB( 21) . VUB( 21 ) . G( 100 ) , DF( 21) . A(21,30), 
c WK(IOOOO) 

INTEGER NRA.NCOLA.NRWK.NRIWK. IGRAD.NDV.NCON, ISTRAT. IPRINT. lOPT, 
c lONED. INFO, NGT 
DOUBLE PRECISION OBJ 

SPECIFICATIONS FOR SUBROUTINES 

EXTERNAL SDTAPER 
EXTERNAL DTAPWT 
EXTERNAL DPTGEN 
EXTERNAL DCOMPAR 
EXTERNAL ADS 

SPECIFICATIONS FOR LOCAL VARIABLES 

INTEGER I.NCALLS 
CHARACTER* 12 NAME I L 



99 



100 



BEGIN EXECUTION 

NCALLS=0 

ADS ARRAY DIMENSIONS 

NRA=21 

NC0LA=30 

NRWK=10000 

NRIWK=2000 

ADS PARAMETERS 

(no user provided gradients, 5 constraints) 

(y/ design variables determined in initial design section) 

IGRAD=0 

NC0N=5 

INPUT INITIAL DESIGN. OUTPUT FILE, 
BOUND ON DESIGN VARIABLES 
OPEN( 88, FILE=' TAPIN.DAT' . STATUS=' OLD' ) 

READ(88,*)NTAP 

NDV=2*(NTAP+1) 

1" < length of segment < 20' 

DO 99 I=1,NTAP 
READ(88.*)TAP(I) 

X(I)=TAP(I) 

VLB(I)=1.0D0/12.0D0 

VUB(I)=20.0D0 

CONTINUE 

0.75" < joint diameter < 10" 

DO 100 I=1,NTAP+1 
READ(88,*)DIAM(I-1) 

X(NTAP+I)=DIAM(I-1) 

VLB(NTAP+I)=0.75D0 

VUB(NTAP+I)=10.0D0 

CONTINUE 

5' < standoff < 30' 

READ (8 8,*) RANGE 
X(2*NTAP+2)=RANGE 
VLB ( 2*NTAP+2 ) =1 5 . ODO 
VUB( 2*NTAP+2 )=30 . ODO 

INPUT LENGTH OF TIME RECORD, TGO 

READ(88,*)TLEN 

READ(88,*)TCO 

INPUT DESIRED PRESSURE PROFILE 

READ( 88 , 1000)NAMFIL 
CLOSE(88) 

OPEN(89,FILE=NAMFIL) 
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non n n n n 



1000 F0RMATCA12) 

0PEN( 88. FILE=' PROFIN.DAT' . STATUS='OLD ' ) 

READ(88.*)NDES 
DO 101 I=0.NDES 

READ(88,*)TDES(I) .PDES(I) 

101 CONTINUE 
CLOSE(88) 

C ID FIVE NONLINEAR CONSTRAINTS 

IDG(1)=0 
IDG(2)=0 
IDG(3)=0 
IDG(A)=0 
IDG(5)=0 

INPUT 

(no ADS print, interractive , optimizer and search options) 
IPRINT=0000 

PRINT*,' ENTER ISTRAT , lOPT , lONED: ' 

READ* , ISTRAT , lOPT , lONED 

OPTIMIZE 

(0 no override. -2 default override) 

INFO=-2 



C 



C 

G 

G 



C 



C 



C 



10 CALL ADS (INFO. ISTRAT. IOPT.IONED,IPRINT.IGRAD.NDV,NCON,X. VLB, 

C VUB , OBJ , G , IDG . NGT . IC , DF , A . NRA , NCOLA , WK , NRWK , IWK , NRIWK) 

IFS TO CONTROL FLOW USING ADS INFO 

IF(INFO.EQ.l) THEN 

EVALUATE OBJECTIVE AND CONSTRAINTS 
ADS VARIABLES TO SDTAPER INPUT 

DO 102 I=1.NTAP 
TAP(I)=X(I) 

102 CONTINUE 

DO 103 I=1,NTAP+1 
DIAM(I-1)=X(NTAP+I) 

103 CONTINUE 
RANGE=X(2*NTAP+2) 

CALCULATE PRESSURE PROFILE 

CALL SDTAPER 
NCALLS=NCALLS+1 

CALCULATE WEIGHT 

CALL DTAPWT 

CORRESPONDING DESIRED PRESSURES 



GALL DPTGEN 

EVALUATE OBJECTIVE FUNCTION 

CALL DCOMPAR 
OBJ=SUMSQ 

PRINT INITIAL DESIGN DATA TO SCREEN 

IF(NCALLS.EO.l) THEN 

PRINT*,' INITIAL DESIGN' 

PRINTlOOl, (TAP(I) ,I=1,NTAP) 

PRINT1002, (DIAM(I) , I=0,NTAP) 

PRINT1003, WEIGHT 
PRINT1005, RANGE 
PRINT1004,OBJ 
END IF 

EVALUATE CONSTRAINTS (25if < W < 125^/) 

G(1)=(25.0D0-WEIGHT) 

G(2)=(WEIGHT-75.0D0) 

EVALUATE CONSTRAINTS ( D1<D2<D3<D4 ) 

G(3)=X(4)-X(5) 

G(4)=X(5)-X(6) 

G(5)=X(6)-X(7) 
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GO TO 10 
ELSE 

IF(INF0.EQ.2) THEN 

C RESERVED FOR ADS GRADIENTS 

GO TO 10 
ELSE 

IFCINFO.EQ. -1) THEN 

C ADS OVERRIDE VALUES 

c (for no scaling, rel FD step, min | FD step|,zero) 

IWK(2)=0 
WK(21)=0.001D0 
WK(22)=0.0001D0 
WK(26)=0.1D0 
WK(37)=1.0D-10 
GO TO 10 
ELSE 

C OPTIMIZATION COMPLETE, INFO=0 

IF(INFO.EQ.O) GO TO 20 
END IF 
END IF 
ENDIF 

20 CONTINUE 

C RECORD RESULTANT PRESSURE PROFILES 

DO 104 I=0.NTIME 

WRITE(89,*) TIME(I) ,PRESS(I) ,PCOMPAR(I) 

104 CONTINUE 
CLOSE(89) 

C PRINT DESIGN DATA TO SCREEN 

PRINT* 

PRINT*,' FINAL DESIGN' 

PRINTlOOl, (TAP(I) , I=1,NTAP) 

1001 FORMAT(' LENGTHS = ',3F10.3,' FT') 

PRINT1002, (DIAM(I) ,I=0,NTAP) 

1002 FORMAT ( ' DIAMETERS = ',4F10.2,' IN') 

PRINT1003, WEIGHT 

1003 FORMAT(' CHARGE WEIGHT =',F6.1,' LB') 

PRINT1005, RANGE 

1005 FORMAT(' RANGE = '.FlO.l,' FT') 

PRINT1004,OBJ 

1004 FORMATC SQRT OF AVG OF (Pdesign-Pdesired) "2 =' , F6.1, 

c' PSI') 

PRINT* 

PRINT1006,NCALLS 

1006 FORMAT(' CALLS TO P-T HISTORY GENERATOR =',I4) 

END 



II. SUBROUTINE DCOMPAR 



c SUBROUTINE DCOMPAR 

c This subroutine, for use with DTAPOPT, computes the square root of 
c of the average sum of the squares of pressure differences. 



c 

c 

c 

c 

c 



PRECISION; 
AUTHOR : 

LAST UPDATE; 
INTERFACE; 



DOUBLE 

William Earl Miller II 
6/17/92 

4 BUSES ; OPTI , OPTDPA 



FORDPTGEN, FORDCOMPAR. 
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c 

c 


VARIABLES 










c 


NAME 


TYPE 


COMMON OR LOCAL 


I/O? 


DESCRIPTION 


c 


I 


I 


local 


na 


local index 


c 


NTIME 


I 


OPTI 


I 


calculated t 


c 


PCOMPAR( 1,1=0. NTIME ) 


DP 


FORDPTGEN 


I 


intp'd desired P 


c 


PRESS (I, 1=0, NTIME) 


DP 


OPTDPA 


I 


calculated P 


c 


SUMSQ 


DP 


FORDCOMPAR 


0 


value to be opt'd 



SUBROUTINE DCOMPAR 
IMPLICIT NONE 



C SPECIFICATIONS FOR GLOBAL VARIABLES 

COMMON/OPT I /NTAP , NT I ME 
INTEGER NTAP,NTIME 

COMMON/OPTDPA/ TAP, DIAM, TIME. PRESS 

DOUBLE PRECISION TAP( 10) , DIAM( 0 : 10 ), TIME ( 0 : 1000) , PRESS ( 0 : 1000 ) 
COMMON/FORDPTGEN/NDES , TDES , PDES , PCOMPAR 

INTEGER NDES 

DOUBLE PRECISION TDES(0:1000) , PDES (0 : 1000) . PCOMPAR(0 : 1000) 

COMMON/FORDCOMPAR/SUMSO 

DOUBLE PRECISION SUMSQ 

C SPECIFICATIONS FOR LOCAL VARIABLE 

INTEGER I 

C BEGIN EXECUTION 

SUMSQ=0 . ODO 
DO 99 I=0,NTIME 

SUMSQ=SUMSO+( PCOMPAR ( I ) - PRESS ( I ) )^-( PCOMPAR ( I ) -PRESS ( I ) ) 

99 CONTINUE 

SUMSQ=SQRT(SUMSOy(NTIME+l)) 

RETURN 

END 



III. SUBROUTINE DPTGEN 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



SUBROUTINE DPTGEN 

This subroutine, for use with DTAPOPT. interpolates the input 
desired pressure -time history to find the desired pressure at each 
time computed by the P-T generator. 



PRECISION: 
AUTHOR : 

LAST UPDATE: 
INTERFACE: 



DOUBLE 

William Earl Miller II 
6/17/92 

3 BUSES: OPTI , OPTDPA, FORDPTGEN 



c VARIABLES 



c 


NAME 


TYPE 


COMMON OR LOCAL 


I/O? 


DESCRIPTION 


c 


I 


I 


local 


na 


local index 


c 


J 


I 


local 


na 


local index 


c 


NDES 


I 


FORDPTGEN 


I 


No. input pts - 1 


c 


NTIME 


I 


OPTI 


I 


No. calc pts - 1 


c 


PCOMPAR(I), 1=0, NTIME 


DP 


FORDPTGEN 


0 


Int Desired Press 


c 


PDES(I), 1=0, NDES 


DP 


FORDPTGEN 


I 


Input pressure 


c 


TDES(I), 1=0, NDES 


DP 


FORDPTGEN 


I 


Input time 


c 


TIME(I), 1=0, NTIME 


DP 


OPTDPA 


I 


Computed time 



SUBROUTINE DPTGEN 
IMPLICIT NONE 

C SPECIFICATIONS FOR GLOBAL VARIABLES 

COMMON/OPTI/NTAP , NTIME 
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INTEGER NTAP.NTIME 

COMMON/OPTDPA/ TAP, DIAM, TIME. PRESS 

DOUBLE PRECISION TAP( 10 ). DIAM( 0 : 10 ), TIME( 0 ; 1000) , PRESS (0 : 1000) 
COMMON/FORDPTGEN/NDES , TDES . PDES . PCOMPAR 

INTEGER NDES 

DOUBLE PRECISION TDES (0 : 1000) , PDES (0 : 1000) , PCOMPAR(0 : 1000) 

C SPECIFICATIONS FOR LOCAL VARIABLES 

INTEGER I,J 

C BEGIN EXECUTION 

J=0 

IF((ABS(TIME(J)) .GE.l.OD-13) .AND. (ABS (TDES ( J) ) . GE . 1 . OD-13 ) ) THEN 
PRINT*, ' INPUT PRESSURE PROFILE START TIME NOT ZERO' 

STOP 

ELSE 

TIME(J)=TDES(J) 

END IF 

DO 100 I=0,NTIME 
1 IF(TIME(I) .GE.TDES(J)) THEN 

IF(TIME(I) .LE.TDES(J+1)) THEN 

PCOMPAR(I)=PDES(J)+((TIME(I)-TDES(J))/(TDES(J+l)-TDES(J))) 
c *(PDES(J+1)-PDES(J)) 

GO TO 100 
ENDIF 
END IF 
J=J+1 

IF(J.EQ.NDES+1) THEN 

PRINT*. 'OUT OF INPUT PLOT POINTS' 

STOP 
ELSE 
GO TO 1 
ENDIF 

100 CONTINUE 
RETURN 
END 



IV. SUBROUTINE DTAPWT 



c 


SUBROUTINE 


DTAPWT 




c 


This 


subroutine , 


for use with DTAPOPT, calculates the weight 


c 

— 


of a tapered charge . 


c 


PRECISION: 


DOUBLE 




c 


AUTHOR : 


William 


Earl Miller II 


c 


LAST UPDATE; 6/17/92 




c 

c - 


INTERFACE: 


3 BUSES: 


: OPTI, OPTDPA, PASS 



c VARIABLES FROM COMMON STATEMENTS 



c 


NAME 


TYPE 


COMMON OR LOCAL 


I/O? 


DESCRIPTION 


c 


DIAM(I) ,I 


=0,NTAP DP 


OPTDPA 


I 


Chg Diams (in) 


c 


I 


I 


local 


na 


local index 


c 


NTAP 


I 


OPTI 


I 


No. Chg Segs 


c 


PI 


DP 


local 


na 


Pi 


c 


RHO 


DP 


local 


na 


H20 dens (lbm/ft'^3) 


c 


SGRAV 


DP 


local 


na 


Sp Grav of Chg 


c 


TAP(I) ,1= 


1 , NTAP DP 


OPTDPA 


I 


Length of Seg (ft) 


c 


VOL 


DP 


local 


na 


Charge Vol 


c 


WEIGHT 


DP 


PASS 


na 


Charge Wt (Ibra) 






SUBROUTINE DTAPWT 
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IMPLICIT NONE 

C SPECIFICATIONS FOR GLOBAL VARIABLES 

COMMON /OPT I /NTAP , NT I ME 
INTEGER NTAP.NTIME 

COMMON/OPTDPA/ TAP, DIAM, TIME, PRESS 

DOUBLE PRECISION TAP( 10 ) , DIAM(0 : 10) , TIME(0 ; 1000 ) , PRESS(0 ; 1000) 
COMMON/PASS/ WEIGHT 
DOUBLE PRECISION WEIGHT 

C SPECIFICATIONS FOR LOCAL VARIABLES 

DOUBLE PRECISION PI ,RHO, SGRAV,VOL 
INTEGER I 

C BEGIN EXECUTION 

PI=4 . 0D0*ATAN( 1 . ODO ) 

RHO= 62.32D0 
SGRAV= 1.712D0 
VOL=0 . ODO 
DO 101 1=1, NTAP 

VOL=VOL+PI*TAP(I)*(DIAM(I-l)*DIAM(I-l)+DIAM(I)*DIAM(I) 
c +DIAM(I-1)*DIAM(I))/1728.0D0 
101 CONTINUE 

WEIGHT=VOL*SGRAV*RHO 

RETURN 

END 
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